/*!
 * \file CFEMNonlinear.h
 *
 * \brief Header file for CFEMNonlinear class
 * \author wanyzh
 * \date Jan 2021
 *
 * 
 */
#pragma once
#include "FEM2DVector.h"
/*!
 * \class CFEMNonlinear
 *
 * \brief FEM class for nonlinear equations
 *
 * \author wanyzh
 * \date Jan 2021
 */
class CFEMNonlinear :
	public CFEM2DVector
{
public:
	CFEMNonlinear(const int &num_vec = 1){
		m_nSizeVec = num_vec;}

	///@brief A general assemble function integration over 2D domain, see also CFEM2DVector::assembleMatrix2D.
	///The different between this version with the version in the base class CFEM2DVector::assembleMatrix2D is that 
	///the coefficient is not a input function which has been defined by the caller.
	///Now the coefficient function is related to the unknowns which means it is a non-linear problem.
	///For general purpose, the design of the coefficient is a function with derivative parameters, that is:
	///\f[
	///\sum_{j=1}^{N_{b}}\int_{E_{n}} \frac{\partial^{d+e} c_{h}}{\partial x^{d} \partial y^{e}} \frac{\partial^{r+s} \psi_{n \alpha}}{\partial x^{r} \partial y^{s}} \frac{\partial^{p+q} \psi_{n \beta}}{\partial x^{p} \partial y^{q}} d x d y
	///\f]
	///The \f$ c_h \f$ in the above integration is the unknown of the element, it is defined by
	///\f[
	///c_h=\sum_{k=1}^{N_{l b}} u_{T_{b}}(k, n) \psi_{n k}
	///\f]
	///And this function has been implemented in the error measurement, see also CFEM2DNonlinear::localFESol.
	///@param uh					The finite coefficients, it is calculated by the FEM and used as a parameters for the next iteration. 
	///@param t						Time for the calculation of the uh.
	///@param basis_uh				The basis test function of the uh.
	///@param basis_der_x_uh		The x derivative of the solution uh.
	///@param basis_der_y_uh		The y derivative of the solution uh.
	///@param Tb					The Tb for the solution uh.
	///@param basis_trial			Basis trial function.		
	///@param basis_trial_der_x		The x derivative of the trial function.
	///@param basis_trial_der_y		The y derivative of the trial function.
	///@param Tb_trial				Tb for trial function.
	///@param basis_test			Basis test function.
	///@param basis_test_der_x		The x derivative of the test function.
	///@param basis_test_der_y		The y derivative of the test function.
	///@param Tb_test				Tb for test function.
	///@param matrix_size			The matrix size get from Pb, it is a pair and the first element represents the rows' number and second one is the cols' number.	
	///Other parameters and return are the same with the function CFEM2D::assembleMatrix2D.
	SpMat assembleMatrix2DNonlinear(Rsv &uh, double t, CBasisFunction &basis_uh, const int &basis_der_x_uh, const int &basis_der_y_uh, vector<CElement>&Tb,
		CBasisFunction &basis_trial, const int &basis_trial_der_x, const int &basis_trial_der_y, vector<CElement>&Tb_trial,
		CBasisFunction &basis_test, const int &basis_test_der_x, const int &basis_test_der_y, vector<CElement>&Tb_test, const pair<int, int> &matrix_size);

	///@brief Local solution over the element.
	///
	///\f[
	///u_h=\sum_{k=1}^{N_{l b}} u_{T_{b}(k, n)} \psi_{n k}
	///\f]
	///@param uh_vector				It is  \f$ u_{T_{b}(k, n)} \f$ in the above equation.
	///@param x						The coordinates of Gauss points.
	///@param vertices				The node coordinates of the element.
	///@param basis_fun				The basis function contains basis type, number of local elements.
	///@param basis_der_x			The x derivative.
	///@param basis_der_y			The y derivative.
	double localFESol(
		const vector<double> &uh_vector,
		const CPoint &x,
		const vector<CPoint> &vertices,
		CBasisFunction &basis_fun,
		const int &basis_der_x,
		const int &basis_der_y);

	///@brief Assemble of vectors for the non-linear term in NS equations.
	///
	///The basic idea is similar to CFEM2DVector::assembleVector, but the load is not a function defined by the caller.
	///It is a term related to the unknowns which means a non-linear problem. 
	///\f[
	///\sum_{j=1}^{N_{b}}\int_{E_{n}} \frac{\partial^{d+e} f_{1 h}}{\partial x^{d} \partial y^{e}} \frac{\partial^{r+s} f_{2 h}}{\partial x^{r} \partial y^{s}} \frac{\partial^{p+q} \psi_{n \beta}}{\partial x^{p} \partial y^{q}} d x d y
	///\f]
	///@param t					Time for the calculation.
	///@param uh1				Finite element function of u1.
	///@param basis_uh1			Basis function for u1.
	///@param basis_der_x_uh1	The x derivative of u1.
	///@param basis_der_y_uh1	The y derivative of u1.
	///@param Tbuh1				The Tb for u1.
	///@param basis_test		The basis test function.
	Rsv assembleVector2DNonlinear(double t,
		Rsv &uh1, CBasisFunction &basis_uh1, const int &basis_der_x_uh1, const int &basis_der_y_uh1, vector<CElement> &Tbuh1, //parameters for the calculation of uh1
		Rsv &uh2, CBasisFunction &basis_uh2, const int &basis_der_x_uh2, const int &basis_der_y_uh2, vector<CElement> &Tbuh2, //parameters for the calculation of uh2
		CBasisFunction & basis_test, const int &basis_der_x_test, const int &basis_der_y_test, vector<CElement> &Tb_test);

protected:
	///@brief Gauss Quadrature over 2D element.
	///
	///This function is  similar with CFEM2DVector::GaussQuad2DCoeTrialTest. The only difference is the coe function's type.
	///\f[
	///\int_{E_{n}} \frac{\partial^{d+e} c_{h}}{\partial x^{d} \partial y^{e}} \frac{\partial^{r+s} \psi_{n \alpha}}{\partial x^{r} \partial y^{s}} \frac{\partial^{p+q} \psi_{n \beta}}{\partial x^{p} \partial y^{q}} d x d y
	///\f]
	///@param basis_trial	The basis trial function.
	///@param basis_test	The basis test function.
	double GaussQuad2DCoeTrialTest(double t, const vector<double>&uh_vector, CBasisFunction &basis_fun_uh, const int &basis_der_x_uh, const int &basis_der_y_uh,
		const vector<double> &Gauss_weights, const vector<CPoint> &Gauss_nodes, const vector<CPoint> &vertices,
		CBasisFunction &basis_trial, const int &basis_index_trial, const int &basis_der_x_trial, const int &basis_der_y_trial,
		CBasisFunction &basis_test, const int &basis_index_test, const int &basis_der_x_test, const int &basis_der_y_test);

	///@brief Gauss quadrature over 2D element.
	///
	///This function is similar with CFEM2DVector::GaussQuad2DLoadTrial. The only difference is the load function's type.
	///\f[
	///\int_{E_{n}} \frac{\partial^{d+e} f_{1 h}}{\partial x^{d} \partial y^{e}} \frac{\partial^{r+s} f_{2 h}}{\partial x^{r} \partial y^{s}} \frac{\partial^{p+q} \psi_{n \beta}}{\partial x^{p} \partial y^{q}} d x d y
	///\f]
	///see also CFEM2DVector::GaussQuad2DLoadTrial and CFEM2DNonlinear::assembleVector2DNonlinear
	double GaussQuad2DLoadTrial(double t,//time for the calculation of uh1 and uh2
		const vector<double> uh1_vector, CBasisFunction &basis_uh1, const int &basis_der_x_uh1, const int &basis_der_y_uh1, //parameters for the calculation of uh1
		const vector<double> uh2_vector, CBasisFunction &basis_uh2, const int &basis_der_x_uh2, const int &basis_der_y_uh2,  //parameters for the calculation of uh2	
		const vector<double> &Gauss_weights, const vector<CPoint> &Gauss_nodes, const vector<CPoint> &vertices,
		CBasisFunction &basis_test, const int &basis_index_test, const int &basis_der_x_test, const int &basis_der_y_test);
};

